Chapter 6 Diversity analysis

6.1 Alpha diversity

# Calculate Hill numbers
richness <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 0) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(richness = 1) %>%
  rownames_to_column(var = "sample")

neutral <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(neutral = 1) %>%
  rownames_to_column(var = "sample")

phylogenetic <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, tree = genome_tree) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(phylogenetic = 1) %>%
  rownames_to_column(var = "sample")

# Aggregate basal GIFT into elements
dist <- genome_gifts %>%
  to.elements(., GIFT_db3) %>%
  traits2dist(., method = "gower")

functional <- genome_counts_filt %>%
  column_to_rownames(var = "genome") %>%
  dplyr::select(where(~ !all(. == 0))) %>%
  hilldiv(., q = 1, dist = dist) %>%
  t() %>%
  as.data.frame() %>%
  dplyr::rename(functional = 1) %>%
  rownames_to_column(var = "sample") %>%
  mutate(functional = if_else(is.nan(functional), 1, functional))

# Merge all metrics
alpha_div <- richness %>%
  full_join(neutral, by = join_by(sample == sample)) %>%
  full_join(phylogenetic, by = join_by(sample == sample)) %>%
  full_join(functional, by = join_by(sample == sample))

6.1.1 Wild samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="0_Wild") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.2 Acclimation samples

alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="1_Acclimation") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = Population, group=Population, color=Population, fill=Population)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#008080', "#d57d2c")) +
      scale_fill_manual(name="Population",
          breaks=c("Cold_wet","Hot_dry"),
          labels=c("Cold","Hot"),
          values=c('#00808050', "#d57d2c50")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.58) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

6.1.3 Post-Transplant_1 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="5_Post-FMT1") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-Cold","Hot-Hot", "Cold-Hot"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
        scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")


# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.1.4 Post-Transplant_2 samples

bxp1<-alpha_div %>%
  pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
  left_join(., sample_metadata, by = join_by(sample == Tube_code)) %>%
  filter(time_point=="6_Post-FMT2") %>%
  mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic","functional"))) %>%
      ggplot(aes(y = value, x = type, group=type, color=type, fill=type)) +
      geom_boxplot(outlier.shape = NA) +
      geom_jitter(alpha=0.5) +
      scale_color_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA","#d57d2c", "#76b183")) +
      scale_fill_manual(name="Type",
          breaks=c("Control","Hot_control", "Treatment"),
          labels=c("Cold-control","Warm-control", "Cold-intervention"),
          values=c("#4477AA50","#d57d2c50","#76b18350")) +
          scale_x_discrete(labels = c("Control" = "Cold-Cold", "Hot_control" = "Hot-Hot", "Treatment" = "Cold-Hot")) +
      facet_wrap(. ~ metric,scales = "free") +
      coord_cartesian(xlim = c(1, NA)) +
      stat_compare_means(size=3, label.x=.7) +
      theme_classic() +
        theme(
    strip.background = element_blank(),
    panel.grid.minor.x = element_line(size = .1, color = "grey"),
    axis.title.x = element_blank(),
    axis.title.y = element_text(size=10),
    axis.text.x = element_text(angle = 45, hjust = 1),
    # Increase plot size
    plot.title = element_text(size = 10),
    axis.text = element_text(size = 8),
    axis.title = element_text(size = 8)
      ) +
  ylab("Alpha diversity")

# Add significance bars with p-values
bxp1 + 
  geom_signif(comparisons = list(c("Control", "Hot_control"), c("Control", "Treatment"), c("Hot_control", "Treatment")),
              map_signif_level = TRUE)

6.2 Beta diversity

beta_q0n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 0)

beta_q1n <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1)

beta_q1p <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, tree = genome_tree)

beta_q1f <- genome_counts_filt %>%
  column_to_rownames(., "genome") %>%
  hillpair(., q = 1, dist = dist)

6.3 Permanovas

6.3.0.1 Load required data

meta <- column_to_rownames(sample_metadata, "Tube_code")

6.3.1 1. Are the wild populations similar?

6.3.1.1 Wild: P.muralis vs P.liolepis

wild <- meta %>%
  filter(time_point == "0_Wild")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

wild.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(wild))]
identical(sort(colnames(wild.counts)), sort(as.character(rownames(wild))))

wild_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild")

6.3.1.1.1 Number of samples used

[1] 27

beta_div_richness_wild<-hillpair(data=wild.counts, q=0)
beta_div_neutral_wild<-hillpair(data=wild.counts, q=1)
beta_div_phylo_wild<-hillpair(data=wild.counts, q=1, tree=genome_tree)
beta_div_func_wild<-hillpair(data=wild.counts, q=1, dist=dist)

6.3.1.2 Richness

betadisper(beta_div_richness_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000012 0.000012 0.0012    999  0.978
Residuals 25 0.257281 0.010291                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.98
Hot_dry   0.97302        

adonis2(formula=beta_div_richness_wild$S ~ Population, data=wild[labels(beta_div_richness_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	1.542719	0.2095041	6.625717	0.001
Residual	25	5.820951	0.7904959	NA	NA
Total	26	7.363669	1.0000000	NA	NA

6.3.1.3 Neutral

betadisper(beta_div_neutral_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.000048 0.0000476 0.0044    999   0.95
Residuals 25 0.270114 0.0108046                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.953
Hot_dry   0.94763        

adonis2(formula=beta_div_neutral_wild$S ~ Population, data=wild[labels(beta_div_neutral_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	1.918266	0.2608511	8.822682	0.001
Residual	25	5.435610	0.7391489	NA	NA
Total	26	7.353876	1.0000000	NA	NA

6.3.1.4 Phylogenetic

betadisper(beta_div_phylo_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03585 0.035847 2.4912    999  0.136
Residuals 25 0.35973 0.014389                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.131
Hot_dry   0.12705        

adonis2(formula=beta_div_phylo_wild$S ~ Population, data=wild[labels(beta_div_phylo_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	0.3218613	0.2162815	6.899207	0.001
Residual	25	1.1662981	0.7837185	NA	NA
Total	26	1.4881594	1.0000000	NA	NA

6.3.1.5 Functional

betadisper(beta_div_func_wild$S, wild$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq     F N.Perm Pr(>F)
Groups     1 0.019387 0.019387 1.653    999  0.199
Residuals 25 0.293200 0.011728                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.205
Hot_dry   0.21033        

adonis2(formula=beta_div_func_wild$S ~ Population, data=wild[labels(beta_div_func_wild$S),], permutations=999) %>%
            as.matrix() %>%
            kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	0.0831048	0.1680538	5.05002	0.067
Residual	25	0.4114083	0.8319462	NA	NA
Total	26	0.4945131	1.0000000	NA	NA

beta_q0n_nmds_wild <- beta_div_richness_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_wild <- beta_div_neutral_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_wild <- beta_div_phylo_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_wild <- beta_div_func_wild$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                vegan::scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(wild_nmds, by = join_by(sample == Tube_code))

6.3.2 2. Effect of acclimation

accli <- meta %>%
  filter(time_point == "1_Acclimation")

# Create a temporary modified version of genome_counts_filt
temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli.counts <- temp_genome_counts[, which(colnames(temp_genome_counts) %in% rownames(accli))]
identical(sort(colnames(accli.counts)), sort(as.character(rownames(accli))))

accli_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation")

6.3.2.0.1 Number of samples used

[1] 26

beta_div_richness_accli<-hillpair(data=accli.counts, q=0)
beta_div_neutral_accli<-hillpair(data=accli.counts, q=1)
beta_div_phylo_accli<-hillpair(data=accli.counts, q=1, tree=genome_tree)
beta_div_func_accli<-hillpair(data=accli.counts, q=1, dist=dist)

6.3.2.1 Richness

betadisper(beta_div_richness_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)   
Groups     1 0.093187 0.093187 11.812    999  0.003 **
Residuals 24 0.189340 0.007889                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
          Cold_wet Hot_dry
Cold_wet             0.004
Hot_dry  0.0021532        

adonis2(formula=beta_div_richness_accli$S ~ Population, data=accli[labels(beta_div_richness_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	1.639630	0.1929088	5.736415	0.001
Residual	24	6.859879	0.8070912	NA	NA
Total	25	8.499509	1.0000000	NA	NA

6.3.2.2 Neutral

betadisper(beta_div_neutral_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.05603 0.056026 4.1918    999   0.05 *
Residuals 24 0.32077 0.013365                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.048
Hot_dry  0.051717        

adonis2(formula=beta_div_neutral_accli$S ~ Population, data=accli[labels(beta_div_neutral_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	1.961192	0.2493171	7.970889	0.001
Residual	24	5.905063	0.7506829	NA	NA
Total	25	7.866255	1.0000000	NA	NA

6.3.2.3 Phylogenetic

betadisper(beta_div_phylo_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.03637 0.036365 2.3087    999  0.149
Residuals 24 0.37804 0.015752                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet            0.129
Hot_dry   0.14172        

adonis2(formula=beta_div_phylo_accli$S ~ Population, data=accli[labels(beta_div_phylo_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	0.2113395	0.1379515	3.84066	0.015
Residual	24	1.3206449	0.8620485	NA	NA
Total	25	1.5319844	1.0000000	NA	NA

6.3.2.4 Functional

betadisper(beta_div_func_accli$S, accli$Population) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00409 0.004087 0.1789    999  0.681
Residuals 24 0.54821 0.022842                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
         Cold_wet Hot_dry
Cold_wet             0.68
Hot_dry   0.67607        

adonis2(formula=beta_div_func_accli$S ~ Population, data=accli[labels(beta_div_func_accli$S),], permutations=999) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	1	0.0040177	0.0059732	0.1442177	0.63
Residual	24	0.6686040	0.9940268	NA	NA
Total	25	0.6726217	1.0000000	NA	NA

beta_q0n_nmds_accli <- beta_div_richness_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1n_nmds_accli <- beta_div_neutral_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE, trace=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1p_nmds_accli <- beta_div_phylo_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

beta_q1f_nmds_accli <- beta_div_func_accli$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  vegan::scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(accli_nmds, by = join_by(sample == Tube_code))

6.3.3 3. Comparison between Wild and Acclimation

accli1 <- meta  %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

accli1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(accli1))]
identical(sort(colnames(accli1.counts)),sort(as.character(rownames(accli1))))

accli1_nmds <- sample_metadata %>%
  filter(time_point == "0_Wild" | time_point == "1_Acclimation")

6.3.3.0.1 Number of samples used

[1] 53

beta_div_richness_accli1<-hillpair(data=accli1.counts, q=0)
beta_div_neutral_accli1<-hillpair(data=accli1.counts, q=1)
beta_div_phylo_accli1<-hillpair(data=accli1.counts, q=1, tree=genome_tree)
beta_div_func_accli1<-hillpair(data=accli1.counts, q=1, dist=dist)

#Arrange of metadata dataframe
accli1_arrange<-accli1[labels(beta_div_neutral_accli1$S),]
accli1_arrange$Population_time <- interaction(accli1_arrange$Population, accli1_arrange$time_point)

6.3.3.0.2 Richness

betadisper(beta_div_richness_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.02724 0.0272402 3.7857    999  0.057 .
Residuals 51 0.36697 0.0071955                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
                0_Wild 1_Acclimation
0_Wild                         0.065
1_Acclimation 0.057212              

adonis2(formula=beta_div_richness_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_richness_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	3.762933	0.2288365	4.846784	0.001
Residual	49	12.680830	0.7711635	NA	NA
Total	52	16.443763	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_richness_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.6396299	5.736415	0.19290877	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.5273732	1.906862	0.05462715	0.006	0.036	.
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	1.5558412	5.190707	0.17782052	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	1.8259388	8.319131	0.24968031	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.3856333	1.736034	0.09788177	0.009	0.054
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.5427188	6.625717	0.20950408	0.001	0.006	*

6.3.3.0.3 Neutral

betadisper(beta_div_neutral_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.0083 0.0083036 0.9803    999   0.35
Residuals 51 0.4320 0.0084705                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              0_Wild 1_Acclimation
0_Wild                       0.339
1_Acclimation 0.3268              

adonis2(formula=beta_div_neutral_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_neutral_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	4.735485	0.2945657	6.820252	0.001
Residual	49	11.340673	0.7054343	NA	NA
Total	52	16.076158	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_neutral_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	1.9611917	7.970889	0.24931708	0.001	0.006	*
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.7466371	3.064914	0.08498327	0.001	0.006	*
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	2.1421077	8.176866	0.25412251	0.001	0.006	*
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	2.0371666	10.078295	0.28730857	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.6314393	3.060027	0.16054681	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	1.9182663	8.822682	0.26085105	0.001	0.006	*

6.3.3.0.4 Phylogenetic

betadisper(beta_div_phylo_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00026 0.0002572 0.0158    999  0.903
Residuals 51 0.82986 0.0162718                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.904
1_Acclimation 0.90044              

adonis2(formula=beta_div_phylo_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_phylo_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	0.7622566	0.2345983	5.006223	0.001
Residual	49	2.4869430	0.7654017	NA	NA
Total	52	3.2491997	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_phylo_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.2113395	3.840660	0.1379515	0.009	0.054
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.2076830	4.216762	0.1133028	0.003	0.018	.
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.3635093	4.799786	0.1666605	0.002	0.012	.
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.2121479	7.924059	0.2406769	0.001	0.006	*
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.2092433	3.885515	0.1953942	0.001	0.006	*
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.3218613	6.899207	0.2162815	0.001	0.006	*

6.3.3.0.5 Functional

betadisper(beta_div_func_accli1$S, accli1$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00006 0.0000604 0.0032    999  0.954
Residuals 51 0.95116 0.0186502                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               0_Wild 1_Acclimation
0_Wild                        0.954
1_Acclimation 0.95485              

adonis2(formula=beta_div_func_accli1$S ~ time_point*Population, data=accli1[labels(beta_div_func_accli1$S),], permutations=999, strata=accli1$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	0.1024301	0.0866259	1.549079	0.338
Residual	49	1.0800123	0.9133741	NA	NA
Total	52	1.1824424	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_func_accli1$S,accli1_arrange$Population_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Cold_wet.1_Acclimation vs Hot_dry.1_Acclimation	1	0.0040176901	0.144217742	0.0059731793	0.646	1.000
Cold_wet.1_Acclimation vs Cold_wet.0_Wild	1	0.0002775823	0.009969841	0.0003020251	0.698	1.000
Cold_wet.1_Acclimation vs Hot_dry.0_Wild	1	0.0871206378	3.762920746	0.1355376396	0.067	0.402
Hot_dry.1_Acclimation vs Cold_wet.0_Wild	1	0.0025235151	0.120315257	0.0047895600	0.595	1.000
Hot_dry.1_Acclimation vs Hot_dry.0_Wild	1	0.0409732757	4.066330477	0.2026444487	0.127	0.762
Cold_wet.0_Wild vs Hot_dry.0_Wild	1	0.0831047968	5.050019518	0.1680537849	0.048	0.288

beta_richness_nmds_accli1 <- beta_div_richness_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_accli1 <- beta_div_neutral_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_accli1 <- beta_div_phylo_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_accli1 <- beta_div_func_accli1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(accli1_nmds, by = join_by(sample == Tube_code))

6.3.4 4. Effect of FMT on microbiota community

6.3.4.1 Comparison between Acclimation vs Post-FMT1

post3 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post3.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post3))]
identical(sort(colnames(post3.counts)),sort(as.character(rownames(post3))))

post3_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "5_Post-FMT1")

6.3.4.1.1 Number of samples used

[1] 52

beta_div_richness_post3<-hillpair(data=post3.counts, q=0)
beta_div_neutral_post3<-hillpair(data=post3.counts, q=1)
beta_div_phylo_post3<-hillpair(data=post3.counts, q=1, tree=genome_tree)
beta_div_func_post3<-hillpair(data=post3.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post3_arrange<-post3[labels(beta_div_neutral_post3$S),]
post3_arrange$type_time <- interaction(post3_arrange$type, post3_arrange$time_point)

6.3.4.1.2 Richness

betadisper(beta_div_richness_post3$S, post3_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)    
Groups     2 0.087391 0.043695 9.5886    999  0.001 ***
Residuals 49 0.223293 0.004557                         
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
               Control Hot_control Treatment
Control                 0.00400000     0.524
Hot_control 0.00219320                 0.001
Treatment   0.49723207  0.00010148          

adonis2(formula=beta_div_richness_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	4.399457	0.2450497	2.986233	0.001
Residual	46	13.553865	0.7549503	NA	NA
Total	51	17.953321	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_richness_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.224	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.002	0.030	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.010	0.150
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.106	1.000

6.3.4.1.3 Neutral

betadisper(beta_div_neutral_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00430 0.0043011 0.5537    999  0.473
Residuals 50 0.38839 0.0077677                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                      0.47
5_Post-FMT1         0.46029            

adonis2(formula=beta_div_neutral_post3$S ~ time_point*type, data=post3[labels(beta_div_neutral_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	5.32464	0.3149916	4.230493	0.001
Residual	46	11.57943	0.6850084	NA	NA
Total	51	16.90407	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_neutral_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.460	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.015	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.001	0.015	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.014	0.210
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.040	0.600

6.3.4.1.4 Phylogenetic

betadisper(beta_div_phylo_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02481 0.024807 1.6411    999  0.236
Residuals 50 0.75581 0.015116                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.226
5_Post-FMT1         0.20609            

adonis2(formula=beta_div_phylo_post3$S ~ time_point*type, data=post3[labels(beta_div_phylo_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	0.7221892	0.2434166	2.959929	0.002
Residual	46	2.2446957	0.7565834	NA	NA
Total	51	2.9668849	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_phylo_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.540	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.008	0.120
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.020	0.300
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.176	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.005	0.075
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.004	0.060
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.003	0.045	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.045	0.675
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.015	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.013	0.195
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.005	0.075
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.784	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.108	1.000
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.694	1.000

6.3.4.1.5 Functional

betadisper(beta_div_func_post3$S, post3$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     1 0.02408 0.024081 1.2619    999  0.251
Residuals 50 0.95415 0.019083                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1
1_Acclimation                     0.245
5_Post-FMT1         0.26666            

adonis2(formula=beta_div_func_post3$S ~ time_point*type, data=post3[labels(beta_div_func_post3$S),], permutations=999,strata=post3$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	0.1702325	0.1167018	1.215508	0.245
Residual	46	1.2884648	0.8832982	NA	NA
Total	51	1.4586973	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_func_post3$S, post3_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.022076664	0.65002192	0.041534889	0.430	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.016749605	0.52091146	0.033561912	0.479	1
Control.1_Acclimation vs Control.5_Post-FMT1	1	-0.008325555	-0.22800110	-0.015434681	0.900	1
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.126527918	3.45550519	0.197960767	0.076	1
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.051008429	1.23951838	0.076327287	0.283	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.001022590	0.05430389	0.003382513	0.662	1
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.002157067	0.09411569	0.005847832	0.618	1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.056602363	2.56037069	0.145803909	0.179	1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.009569124	0.35095521	0.021463896	0.495	1
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.001745663	-0.08225018	-0.005167199	0.703	1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.057758674	2.84545622	0.159449901	0.148	1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.005575266	0.21803560	0.013444020	0.559	1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.119540855	4.84764704	0.244242909	0.079	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.052587837	1.77308932	0.099762584	0.215	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.012980354	0.44307662	0.028690955	0.466	1

beta_richness_nmds_post3 <- beta_div_richness_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post3 <- beta_div_neutral_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post3 <- beta_div_phylo_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post3 <- beta_div_func_post3$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post3_nmds, by = join_by(sample == Tube_code))

6.3.4.2 Comparison between Acclimation vs Post-FMT2

post4 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post4.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post4))]
identical(sort(colnames(post4.counts)),sort(as.character(rownames(post4))))

post4_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2")

6.3.4.2.1 Number of samples used

[1] 53

beta_div_richness_post4<-hillpair(data=post4.counts, q=0)
beta_div_neutral_post4<-hillpair(data=post4.counts, q=1)
beta_div_phylo_post4<-hillpair(data=post4.counts, q=1, tree=genome_tree)
beta_div_func_post4<-hillpair(data=post4.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post4_arrange<-post4[labels(beta_div_neutral_post4$S),]
post4_arrange$type_time <- interaction(post4_arrange$type, post4_arrange$time_point)

6.3.4.2.2 Richness

betadisper(beta_div_richness_post4$S, post4_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)  
Groups     2 0.05757 0.0287844 3.754    999  0.031 *
Residuals 50 0.38338 0.0076677                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0840000     0.497
Hot_control 0.0769856                 0.003
Treatment   0.4949731   0.0047188          

adonis2(formula=beta_div_richness_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	3.279021	0.1937289	3.924535	0.001
Residual	49	13.646799	0.8062711	NA	NA
Total	52	16.925820	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_richness_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.251	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.027	0.405
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.004	0.060
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.002	0.030	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.034	0.510
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.4.2.3 Neutral

betadisper(beta_div_neutral_post4$S, post4_arrange$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00848 0.0084800 1.2368    999  0.269
Residuals 51 0.34967 0.0068563                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                      0.26
6_Post-FMT2          0.2713            

adonis2(formula=beta_div_neutral_post4$S ~ time_point*Population, data=post4[labels(beta_div_neutral_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	3.823089	0.2357666	5.038847	0.001
Residual	49	12.392477	0.7642334	NA	NA
Total	52	16.215567	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_neutral_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.479	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.015	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.015	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.006	0.090
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.015	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.002	0.030	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.015	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.030	0.450
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.015	.

6.3.4.2.4 Phylogenetic

betadisper(beta_div_phylo_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     1 0.03769 0.037692 3.8614    999  0.045 *
Residuals 51 0.49783 0.009761                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                      0.04
6_Post-FMT2        0.054869            

adonis2(formula=beta_div_phylo_post4$S ~ time_point*Population, data=post4[labels(beta_div_phylo_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	0.6466392	0.2442162	5.277784	0.001
Residual	49	2.0011759	0.7557838	NA	NA
Total	52	2.6478151	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_phylo_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.789330	0.04999136	0.548	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.054731	0.16919285	0.015	0.225
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.379209	0.22597458	0.005	0.075
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.879364	0.11134092	0.124	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.878754	0.24542555	0.003	0.045	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.052184	0.20208192	0.001	0.015	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.396667	0.28560797	0.001	0.015	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.829222	0.19311005	0.015	0.225
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.846334	0.26761166	0.001	0.015	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.383624	0.21505615	0.002	0.030	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.060129	0.24027055	0.001	0.015	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.499774	0.13512455	0.021	0.315
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.382025	0.12958449	0.028	0.420
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.722460	0.14541146	0.007	0.105
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.043656	0.20174244	0.001	0.015	.

6.3.4.2.5 Functional

betadisper(beta_div_func_post4$S, post4$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     1 0.00024 0.0002438 0.0119    999  0.925
Residuals 51 1.04336 0.0204581                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 6_Post-FMT2
1_Acclimation                     0.926
6_Post-FMT2         0.91349            

adonis2(formula=beta_div_func_post4$S ~ time_point*Population, data=post4[labels(beta_div_func_post4$S),], permutations=999,strata=post4$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	3	0.0258537	0.0188483	0.313769	0.578
Residual	49	1.3458235	0.9811517	NA	NA
Total	52	1.3716772	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_func_post4$S, post4_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.0220766636	0.650021923	0.0415348890	0.407	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.0167496054	0.520911459	0.0335619117	0.478	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0376809680	1.119650309	0.0694587220	0.307	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.0292200956	0.920511083	0.0578191917	0.390	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0133433458	0.270473784	0.0177122064	0.510	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.0010225901	0.054303886	0.0033825127	0.630	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-0.0005177706	-0.025585400	-0.0016016487	0.727	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.0013301207	0.072110871	0.0044867082	0.600	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0060959077	0.174487757	0.0107878382	0.613	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0010345754	0.055797964	0.0034752533	0.628	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-0.0001056284	-0.006306177	-0.0003942915	0.689	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0017235602	0.051851181	0.0032302306	0.755	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.442986255	-0.0284750185	0.851	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.034047378	-0.0021324990	0.904	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.110487573	0.0068581148	0.712	1

beta_richness_nmds_post4 <- beta_div_richness_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post4 <- beta_div_neutral_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post4 <- beta_div_phylo_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post4 <- beta_div_func_post4$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post4_nmds, by = join_by(sample == Tube_code))

6.3.4.3 Comparison between Acclimation vs Post-FMT1 and Post-FMT2

post6 <- meta  %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

temp_genome_counts <- transform(genome_counts_filt, row.names = genome_counts_filt$genome)
temp_genome_counts$genome <- NULL

post6.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post6))]
identical(sort(colnames(post6.counts)),sort(as.character(rownames(post6))))

post6_nmds <- sample_metadata %>%
  filter(time_point == "1_Acclimation" | time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

6.3.4.3.1 Number of samples used

[1] 79

beta_div_richness_post6<-hillpair(data=post6.counts, q=0)
beta_div_neutral_post6<-hillpair(data=post6.counts, q=1)
beta_div_phylo_post6<-hillpair(data=post6.counts, q=1, tree=genome_tree)
beta_div_func_post6<-hillpair(data=post6.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post6_arrange<-post6[labels(beta_div_neutral_post6$S),]
post6_arrange$type_time <- interaction(post6_arrange$type, post6_arrange$time_point)

6.3.4.3.2 Richness

betadisper(beta_div_richness_post6$S, post6_arrange$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.06454 0.032268 4.9753    999  0.014 *
Residuals 76 0.49290 0.006485                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0120000     0.681
Hot_control 0.0119513                 0.003
Treatment   0.6817089   0.0035106          

adonis2(formula=beta_div_richness_post6$S ~ time_point*type, data=post6[labels(beta_div_richness_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	8	6.395467	0.2451322	2.841433	0.001
Residual	70	19.694403	0.7548678	NA	NA
Total	78	26.089870	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_richness_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.3657243	1.123239	0.06966584	0.244	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.2605860	4.943717	0.24788342	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.6804988	2.100611	0.12283837	0.001	0.036	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8481258	2.673290	0.16033371	0.001	0.036	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.1180661	3.809275	0.20252111	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.8072604	2.940901	0.16392161	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4817387	1.660775	0.09968176	0.029	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.1179704	3.885459	0.20573812	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.3606630	5.087152	0.24124415	0.001	0.036	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.7216200	2.172734	0.11956009	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9551308	2.926054	0.16322910	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.2263345	4.039487	0.20157637	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.9130048	3.195028	0.16645080	0.002	0.072
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.5959230	1.984036	0.11032208	0.002	0.072
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.2747787	4.275366	0.21086503	0.002	0.072
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.4319792	5.384836	0.25180628	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.8172413	3.194690	0.17558364	0.002	0.072
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.5796135	2.441615	0.13239702	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6397330	2.913695	0.15405213	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.4575447	6.224524	0.28007456	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3276169	1.412318	0.08111028	0.036	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.011	0.396
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.036	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.001	0.036	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.266	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.110	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.122	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.008	0.288
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.002	0.072
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.466	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.002	0.072
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.036	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.036	.

6.3.4.3.3 Neutral

betadisper(beta_div_neutral_post6$S, post6_arrange$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00896 0.0044799 0.6329    999  0.508
Residuals 76 0.53799 0.0070788                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                   0.45700       0.251
5_Post-FMT1         0.45960                   0.756
6_Post-FMT2         0.27047     0.75275            

adonis2(formula=beta_div_neutral_post6$S ~ time_point*type, data=post6[labels(beta_div_neutral_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	8	7.736408	0.3121974	3.971673	0.001
Residual	70	17.044094	0.6878026	NA	NA
Total	78	24.780501	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_neutral_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.2759858	0.9928976	0.06208366	0.453	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	1.4146464	6.5078325	0.30257965	0.001	0.036	.
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.8153894	3.0970603	0.17113610	0.001	0.036	.
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	1.1809241	4.4856470	0.24265567	0.001	0.036	.
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.4321524	5.7774260	0.27806264	0.001	0.036	.
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	1.1524353	4.9536068	0.24825621	0.001	0.036	.
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.4748999	1.9609749	0.11561688	0.004	0.144
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.3503168	5.4081420	0.26499923	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	1.6347704	6.8326887	0.29925029	0.001	0.036	.
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.9517634	3.3715700	0.17404733	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	1.3127773	4.6298256	0.23585668	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	1.6713369	6.2395460	0.28056085	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	1.3540292	5.3398081	0.25022756	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.6311089	2.4041625	0.13063146	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	1.6125755	5.9825981	0.27215155	0.001	0.036	.
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	1.5409781	6.8338056	0.29928456	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.9133614	4.0964534	0.21451383	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.6954835	3.2951234	0.17077493	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.6202327	3.1519868	0.16457754	0.001	0.036	.
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	1.5701179	7.6327037	0.32297209	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.3634438	1.7083388	0.09647087	0.041	1.000
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.2508491	0.13047758	0.011	0.396
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.1436369	0.20570451	0.001	0.036	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.7146920	0.18842252	0.002	0.072
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.5521758	0.08843210	0.085	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.1302726	0.24279254	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.6372683	0.09840968	0.045	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.3252026	0.08117526	0.152	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.2700097	0.22158835	0.001	0.036	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.9200493	0.16294873	0.002	0.072
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.9421261	0.10824392	0.009	0.324
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.8752789	0.26858075	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.3161367	0.07600637	0.169	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.6483346	0.22511910	0.001	0.036	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.2065321	0.12119453	0.003	0.108
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.7710313	0.26507845	0.001	0.036	.

6.3.4.3.4 Phylogenetic

betadisper(beta_div_phylo_post6$S, post6$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.04250 0.021252 1.8628    999  0.161
Residuals 76 0.86703 0.011408                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                  0.223000       0.051
5_Post-FMT1        0.205952                   0.746
6_Post-FMT2        0.054864    0.717296            

adonis2(formula=beta_div_phylo_post6$S ~ time_point*type, data=post6[labels(beta_div_phylo_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	8	1.055866	0.2692287	3.223651	0.001
Residual	70	2.865952	0.7307713	NA	NA
Total	78	3.921818	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_phylo_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.05690611	0.7893300	0.04999136	0.545	1.000
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.10866997	3.0547314	0.16919285	0.010	0.360
Control.1_Acclimation vs Control.5_Post-FMT1	1	0.11760287	2.9988032	0.16661126	0.024	0.864
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09624422	1.7247115	0.10968160	0.163	1.000
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.18586788	4.1904227	0.21836011	0.003	0.108
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.15591805	4.3792085	0.22597458	0.005	0.180
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.07099742	1.8793639	0.11134092	0.120	1.000
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.18682367	4.8787543	0.24542555	0.001	0.036	.
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.23108846	4.0521838	0.20208192	0.003	0.108
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.26358465	4.3608960	0.21417997	0.002	0.072
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.25319427	3.2738422	0.17915456	0.044	1.000
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.39050120	5.9837393	0.27218933	0.001	0.036	.
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	0.36496892	6.3966666	0.28560797	0.001	0.036	.
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.22628210	3.8292220	0.19311005	0.015	0.540
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.34830814	5.8463335	0.26761166	0.002	0.072
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	0.14203376	5.4200212	0.25303529	0.001	0.036	.
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.09666753	2.3682173	0.13635351	0.014	0.504
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.09252600	2.9824958	0.15711821	0.012	0.432
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.10002871	4.3836237	0.21505615	0.002	0.072
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	0.12577510	5.0601287	0.24027055	0.001	0.036	.
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.06334378	2.4997737	0.13512455	0.014	0.504
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.795	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.102	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.006	0.216
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.162	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.002	0.072
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.704	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.132	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.295	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.019	0.684
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.035	1.000
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.036	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.089	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.028	1.000
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.004	0.144
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.002	0.072

6.3.4.3.5 Functional

betadisper(beta_div_func_post6$S, post6$time_point) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq  Mean Sq      F N.Perm Pr(>F)
Groups     2 0.03015 0.015074 0.7893    999  0.487
Residuals 76 1.45147 0.019098                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              1_Acclimation 5_Post-FMT1 6_Post-FMT2
1_Acclimation                   0.29300       0.932
5_Post-FMT1         0.26551                   0.298
6_Post-FMT2         0.93087     0.28011            

adonis2(formula=beta_div_func_post6$S ~ time_point*type, data=post6[labels(beta_div_func_post6$S),], permutations=999,strata=post6$individual) %>%
  as.matrix() %>%
  kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	8	0.1632318	0.0763851	0.7236452	0.472
Residual	70	1.9737271	0.9236149	NA	NA
Total	78	2.1369589	1.0000000	NA	NA

pairwise <- pairwise.adonis(beta_div_func_post6$S, post6_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.1_Acclimation vs Treatment.1_Acclimation	1	0.0220766636	0.650021923	0.0415348890	0.456	1
Control.1_Acclimation vs Hot_control.1_Acclimation	1	0.0167496054	0.520911459	0.0335619117	0.483	1
Control.1_Acclimation vs Control.5_Post-FMT1	1	-0.0083255553	-0.228001097	-0.0154346814	0.882	1
Control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.1265279179	3.455505193	0.1979607668	0.077	1
Control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0510084290	1.239518381	0.0763272870	0.297	1
Control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0376809680	1.119650309	0.0694587220	0.304	1
Control.1_Acclimation vs Control.6_Post-FMT2	1	0.0292200956	0.920511083	0.0578191917	0.374	1
Control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0133433458	0.270473784	0.0177122064	0.505	1
Treatment.1_Acclimation vs Hot_control.1_Acclimation	1	0.0010225901	0.054303886	0.0033825127	0.645	1
Treatment.1_Acclimation vs Control.5_Post-FMT1	1	0.0021570675	0.094115687	0.0058478321	0.635	1
Treatment.1_Acclimation vs Treatment.5_Post-FMT1	1	0.0566023632	2.560370692	0.1458039091	0.161	1
Treatment.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0095691236	0.350955208	0.0214638964	0.509	1
Treatment.1_Acclimation vs Treatment.6_Post-FMT2	1	-0.0005177706	-0.025585400	-0.0016016487	0.715	1
Treatment.1_Acclimation vs Control.6_Post-FMT2	1	0.0013301207	0.072110871	0.0044867082	0.592	1
Treatment.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0060959077	0.174487757	0.0107878382	0.596	1
Hot_control.1_Acclimation vs Control.5_Post-FMT1	1	-0.0017456629	-0.082250179	-0.0051671989	0.719	1
Hot_control.1_Acclimation vs Treatment.5_Post-FMT1	1	0.0577586745	2.845456220	0.1594499006	0.146	1
Hot_control.1_Acclimation vs Hot_control.5_Post-FMT1	1	0.0055752661	0.218035597	0.0134440201	0.548	1
Hot_control.1_Acclimation vs Treatment.6_Post-FMT2	1	0.0010345754	0.055797964	0.0034752533	0.659	1
Hot_control.1_Acclimation vs Control.6_Post-FMT2	1	-0.0001056284	-0.006306177	-0.0003942915	0.673	1
Hot_control.1_Acclimation vs Hot_control.6_Post-FMT2	1	0.0017235602	0.051851181	0.0032302306	0.741	1
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.1195408549	4.847647043	0.2442429086	0.072	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0525878365	1.773089316	0.0997625840	0.226	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0265995825	1.175418056	0.0684360667	0.297	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0145818992	0.699759916	0.0419023938	0.414	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0080695208	-0.216173226	-0.0136958691	0.931	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0129803540	0.443076619	0.0286909552	0.461	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0267162134	1.225605811	0.0755352882	0.272	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0384388433	1.932815819	0.1141461550	0.241	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.0553988290	1.478193905	0.0897060633	0.230	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-0.0040061386	-0.148504693	-0.0093684974	0.744	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0024023972	0.095389804	0.0059265296	0.629	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0004960759	-0.011903277	-0.0007445087	0.849	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.442986255	-0.0284750185	0.860	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.034047378	-0.0021324990	0.917	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.110487573	0.0068581148	0.686	1

beta_richness_nmds_post6 <- beta_div_richness_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = c("sample" = "Tube_code"))

beta_neutral_nmds_post6 <- beta_div_neutral_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = c("sample" = "Tube_code"))

beta_phylo_nmds_post6 <- beta_div_phylo_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = join_by(sample == Tube_code))

beta_func_nmds_post6 <- beta_div_func_post6$S %>%
  metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
  scores() %>%
  as_tibble(., rownames = "sample") %>%
  left_join(post6_nmds, by = join_by(sample == Tube_code))

6.3.5 5. Are there differences between the control and the treatment group?

6.3.5.1 After 1 week –> Post-FMT1

post1 <- meta %>%
  filter(time_point == "5_Post-FMT1")

post1.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post1))]
identical(sort(colnames(post1.counts)),sort(as.character(rownames(post1))))

post1_nmds <- sample_metadata %>%
  filter(time_point == "5_Post-FMT1")

6.3.5.1.1 Number of samples used

[1] 26

beta_div_richness_post1<-hillpair(data=post1.counts, q=0)
beta_div_neutral_post1<-hillpair(data=post1.counts, q=1)
beta_div_phylo_post1<-hillpair(data=post1.counts, q=1, tree=genome_tree)
beta_div_func_post1<-hillpair(data=post1.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post1_arrange<-post1[labels(beta_div_neutral_post1$S),]

6.3.5.1.2 Richness

betadisper(beta_div_richness_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)  
Groups     2 0.017675 0.0088373 2.3825    999  0.097 .
Residuals 23 0.085312 0.0037092                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
              Control Hot_control Treatment
Control                 0.0040000     0.647
Hot_control 0.0068795                 0.231
Treatment   0.6248469   0.2084296          

adonis2(formula=beta_div_richness_post1$S ~ type, data=post1[labels(beta_div_richness_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	1.195567	0.1448246	1.947534	0.001
Residual	23	7.059710	0.8551754	NA	NA
Total	25	8.255277	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_richness_post1$S,post1_arrange$type, perm=999)
pairwise

                     pairs Df SumsOfSqs  F.Model        R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.5615418 1.729004 0.1033537   0.010      0.030   .
2   Control vs Hot_control  1 0.8438429 2.793772 0.1486541   0.001      0.003   *
3 Treatment vs Hot_control  1 0.3734921 1.268929 0.0779971   0.108      0.324    

6.3.5.1.3 Neutral

betadisper(beta_div_neutral_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.011001 0.0055005 0.6303    999   0.58
Residuals 23 0.200714 0.0087267                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.20200     0.957
Hot_control 0.21166                 0.445
Treatment   0.95468     0.43604          

adonis2(formula=beta_div_neutral_post1$S ~ type, data=post1[labels(beta_div_neutral_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	1.395968	0.1900228	2.697931	0.001
Residual	23	5.950350	0.8099772	NA	NA
Total	25	7.346318	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_neutral_post1$S,post1_arrange$type, perm=999)
pairwise

                     pairs Df SumsOfSqs  F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.6051778 2.250849 0.13047758   0.012      0.036   .
2   Control vs Hot_control  1 1.0528902 4.143637 0.20570451   0.001      0.003   *
3 Treatment vs Hot_control  1 0.4150076 1.637268 0.09840968   0.061      0.183    

6.3.5.1.4 Phylogenetic

betadisper(beta_div_phylo_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00440 0.0021994 0.1369    999  0.911
Residuals 23 0.36941 0.0160614                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.92200     0.692
Hot_control 0.91505                 0.770
Treatment   0.63312     0.73046          

adonis2(formula=beta_div_phylo_post1$S ~ type, data=post1[labels(beta_div_phylo_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	0.0745104	0.0705947	0.8735033	0.537
Residual	23	0.9809570	0.9294053	NA	NA
Total	25	1.0554673	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_phylo_post1$S,post1_arrange$type, perm=999)
pairwise

                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.01842535 0.4144162 0.02688498   0.771       1.00    
2   Control vs Hot_control  1 0.05987967 1.7387847 0.09802164   0.110       0.33    
3 Treatment vs Hot_control  1 0.03212966 0.6477782 0.04139746   0.694       1.00    

6.3.5.1.5 Functional

betadisper(beta_div_func_post1$S, post1$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00400 0.0019999 0.1431    999  0.859
Residuals 23 0.32135 0.0139717                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.59500     0.744
Hot_control 0.60188                 0.813
Treatment   0.74597     0.84473          

adonis2(formula=beta_div_func_post1$S ~ type, data=post1[labels(beta_div_func_post1$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	0.1230554	0.1608583	2.204479	0.174
Residual	23	0.6419374	0.8391417	NA	NA
Total	25	0.7649929	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_func_post1$S,post1_arrange$type, perm=999)
pairwise

                     pairs Df  SumsOfSqs   F.Model         R2 p.value p.adjusted sig
1     Control vs Treatment  1 0.11954085 4.8476470 0.24424291   0.068      0.204    
2   Control vs Hot_control  1 0.05258784 1.7730893 0.09976258   0.227      0.681    
3 Treatment vs Hot_control  1 0.01298035 0.4430766 0.02869096   0.453      1.000    

beta_richness_nmds_post1 <- beta_div_richness_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post1 <- beta_div_neutral_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post1 <- beta_div_phylo_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post1 <- beta_div_func_post1$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post1_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post1 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.2 After 2 weeks –>Post-FMT2

post2 <- meta %>%
  filter(time_point == "6_Post-FMT2")

post2.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post2))]
identical(sort(colnames(post2.counts)),sort(as.character(rownames(post2))))

post2_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2")

6.3.5.2.1 Number of samples used

[1] 27

beta_div_richness_post2<-hillpair(data=post2.counts, q=0)
beta_div_neutral_post2<-hillpair(data=post2.counts, q=1)
beta_div_phylo_post2<-hillpair(data=post2.counts, q=1, tree=genome_tree)
beta_div_func_post2<-hillpair(data=post2.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post2_arrange<-post2[labels(beta_div_neutral_post2$S),]

6.3.5.2.2 Richness

betadisper(beta_div_richness_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.002011 0.0010056 0.1982    999  0.845
Residuals 24 0.121775 0.0050740                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.72000     0.818
Hot_control 0.67789                 0.636
Treatment   0.79246     0.59820          

adonis2(formula=beta_div_richness_post2$S ~ type, data=post2[labels(beta_div_richness_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	1.504341	0.1967776	2.939822	0.001
Residual	24	6.140538	0.8032224	NA	NA
Total	26	7.644879	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_richness_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.6463814	2.560441	0.1379515	0.001	0.003	*
Treatment vs Hot_control	1	0.4796256	1.916520	0.1069694	0.001	0.003	*
Control vs Hot_control	1	1.1305044	4.268317	0.2105906	0.001	0.003	*

6.3.5.2.3 Neutral

betadisper(beta_div_neutral_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.008262 0.0041311 0.8024    999  0.434
Residuals 24 0.123559 0.0051483                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.43400     0.651
Hot_control 0.44675                 0.236
Treatment   0.65989     0.25095          

adonis2(formula=beta_div_neutral_post2$S ~ type, data=post2[labels(beta_div_neutral_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	1.923807	0.2603795	4.224537	0.001
Residual	24	5.464666	0.7396205	NA	NA
Total	26	7.388473	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_neutral_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	1.0227481	4.648335	0.2251191	0.001	0.003	*
Treatment vs Hot_control	1	0.5010202	2.206532	0.1211945	0.001	0.003	*
Control vs Hot_control	1	1.3619424	5.771031	0.2650785	0.001	0.003	*

6.3.5.2.4 Phylogenetic

betadisper(beta_div_phylo_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.000407 0.0002034 0.0487    999  0.956
Residuals 24 0.100305 0.0041794                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.93900     0.834
Hot_control 0.93765                 0.771
Treatment   0.83933     0.76015          

adonis2(formula=beta_div_phylo_post2$S ~ type, data=post2[labels(beta_div_phylo_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	0.1594363	0.2042241	3.079623	0.001
Residual	24	0.6212564	0.7957759	NA	NA
Total	26	0.7806927	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_phylo_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	0.05927454	2.382025	0.1295845	0.033	0.099
Treatment vs Hot_control	1	0.06906280	2.722460	0.1454115	0.003	0.009	*
Control vs Hot_control	1	0.11081709	4.043656	0.2017424	0.003	0.009	*

6.3.5.2.5 Functional

betadisper(beta_div_func_post2$S, post2$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01259 0.0062962 0.3249    999  0.777
Residuals 24 0.46507 0.0193778                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.50600     0.646
Hot_control 0.45381                 0.776
Treatment   0.57452     0.74365          

adonis2(formula=beta_div_func_post2$S ~ type, data=post2[labels(beta_div_func_post2$S),], permutations=999) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	2	-0.0037283	-0.0054704	-0.065288	0.923
Residual	24	0.6852623	1.0054704	NA	NA
Total	26	0.6815340	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_func_post2$S,post2_arrange$type, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Treatment vs Control	1	-0.008042888	-0.44298625	-0.028475019	0.847	1
Treatment vs Hot_control	1	-0.001179626	-0.03404738	-0.002132499	0.894	1
Control vs Hot_control	1	0.003630084	0.11048757	0.006858115	0.684	1

beta_richness_nmds_post2 <- beta_div_richness_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post2 <- beta_div_neutral_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post2 <- beta_div_phylo_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post2 <- beta_div_func_post2$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post2_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post2 %>%
            group_by(type) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")

6.3.5.3 Post1 vs Post2

post5 <- meta %>%
  filter(time_point == "6_Post-FMT2" | time_point == "5_Post-FMT1")

post5.counts <- temp_genome_counts[,which(colnames(temp_genome_counts) %in% rownames(post5))]
identical(sort(colnames(post5.counts)),sort(as.character(rownames(post5))))

post5_nmds <- sample_metadata %>%
  filter(time_point == "6_Post-FMT2"| time_point == "5_Post-FMT1")

6.3.5.3.1 Number of samples used

[1] 53

beta_div_richness_post5<-hillpair(data=post5.counts, q=0)
beta_div_neutral_post5<-hillpair(data=post5.counts, q=1)
beta_div_phylo_post5<-hillpair(data=post5.counts, q=1, tree=genome_tree)
beta_div_func_post5<-hillpair(data=post5.counts, q=1, dist=dist)

#Arrange of metadata dataframe
post5_arrange<-post5[labels(beta_div_neutral_post5$S),]
post5_arrange$type_time <- interaction(post5_arrange$type, post5_arrange$time_point)

6.3.5.3.2 Richness

betadisper(beta_div_richness_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.01841 0.0092048 1.7364    999  0.186
Residuals 50 0.26505 0.0053010                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
             Control Hot_control Treatment
Control                 0.041000     0.715
Hot_control 0.039117                 0.228
Treatment   0.716358    0.218648          

adonis2(formula=beta_div_richness_post5$S ~ type*time_point, data=post5[labels(beta_div_richness_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	3.261916	0.1981463	2.322836	0.005
Residual	47	13.200248	0.8018537	NA	NA
Total	52	16.462165	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_richness_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.5615418	1.729004	0.10335366	0.014	0.210
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.8438429	2.793772	0.14865413	0.001	0.015	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.7628135	2.683925	0.14364890	0.001	0.015	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3432605	1.148733	0.06698647	0.270	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.1269580	3.799256	0.19188884	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.3734921	1.268929	0.07799710	0.111	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3571397	1.297184	0.07959561	0.137	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.7769467	2.670898	0.15114670	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.6502360	2.253407	0.13060650	0.003	0.045	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4132091	1.616138	0.09174188	0.003	0.045	.
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0163992	3.760571	0.19030682	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.2732563	1.019281	0.05988979	0.423	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.6463814	2.560441	0.13795154	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.4796256	1.916520	0.10696943	0.001	0.015	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.1305044	4.268317	0.21059061	0.001	0.015	.

6.3.5.3.3 Neutral

betadisper(beta_div_neutral_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq     F N.Perm Pr(>F)
Groups     2 0.01992 0.0099587 1.565    999   0.23
Residuals 50 0.31818 0.0063636                    

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.10800     0.884
Hot_control 0.10701                 0.178
Treatment   0.87156     0.17449          

adonis2(formula=beta_div_neutral_post5$S ~ type*time_point, data=post5[labels(beta_div_neutral_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	3.947979	0.2569798	3.251069	0.002
Residual	47	11.415016	0.7430202	NA	NA
Total	52	15.362995	1.0000000	NA	NA

pairwise<-pairwise.adonis(beta_div_neutral_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.6051778	2.250849	0.13047758	0.020	0.300
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	1.0528902	4.143637	0.20570451	0.002	0.030	.
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.8908158	3.714692	0.18842252	0.001	0.015	.
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.3860927	1.552176	0.08843210	0.097	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	1.3122237	5.130273	0.24279254	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.4150076	1.637268	0.09840968	0.043	0.645
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.3157079	1.325203	0.08117526	0.156	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.0579520	4.270010	0.22158835	0.001	0.015	.
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.7454015	2.920049	0.16294873	0.002	0.030	.
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.4377161	1.942126	0.10824392	0.006	0.090
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	1.3766597	5.875279	0.26858075	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.3176516	1.316137	0.07600637	0.181	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	1.0227481	4.648335	0.22511910	0.001	0.015	.
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.5010202	2.206532	0.12119453	0.002	0.030	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	1.3619424	5.771031	0.26507845	0.001	0.015	.

6.3.5.3.4 Phylogenetic

betadisper(beta_div_phylo_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00051 0.0002543 0.0265    999  0.985
Residuals 50 0.47996 0.0095993                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.89500     0.855
Hot_control 0.88926                 0.930
Treatment   0.82391     0.91902          

adonis2(formula=beta_div_phylo_post5$S ~ type*time_point, data=post5[labels(beta_div_phylo_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	0.3518222	0.180049	2.0641	0.046
Residual	47	1.6022134	0.819951	NA	NA
Total	52	1.9540356	1.000000	NA	NA

pairwise<-pairwise.adonis(beta_div_phylo_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.01842535	0.4144162	0.02688498	0.769	1.000
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.05987967	1.7387847	0.09802164	0.108	1.000
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07917244	3.0180046	0.15869197	0.005	0.075
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.04335491	1.5335604	0.08746429	0.193	1.000
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.10783045	3.7500438	0.18987521	0.002	0.030	.
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.03212966	0.6477782	0.04139746	0.698	1.000
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.06393539	1.5651817	0.09448624	0.173	1.000
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.05265949	1.2240203	0.07544494	0.299	1.000
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.09753501	2.2402429	0.12994265	0.010	0.150
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.07228545	2.3279593	0.12701683	0.037	0.555
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.11759094	3.5538444	0.18174658	0.001	0.015	.
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.06667255	1.9859527	0.11041687	0.101	1.000
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	0.05927454	2.3820253	0.12958449	0.024	0.360
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.06906280	2.7224602	0.14541146	0.003	0.045	.
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.11081709	4.0436561	0.20174244	0.003	0.045	.

6.3.5.3.5 Functional

betadisper(beta_div_func_post5$S, post5$type) %>% permutest(., pairwise = TRUE)

Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999

Response: Distances
          Df  Sum Sq   Mean Sq      F N.Perm Pr(>F)
Groups     2 0.00785 0.0039232 0.2322    999  0.811
Residuals 50 0.84483 0.0168966                     

Pairwise comparisons:
(Observed p-value below diagonal, permuted p-value above diagonal)
            Control Hot_control Treatment
Control                 0.54500     0.623
Hot_control 0.52384                 0.835
Treatment   0.58787     0.85068          

adonis2(formula=beta_div_func_post5$S ~ type*time_point, data=post5[labels(beta_div_func_post5$S),], permutations=999,strata=post5$individual) %>%
        as.matrix() %>%
        kable()

	Df	SumOfSqs	R2	F	Pr(>F)
Model	5	0.1076682	0.075037	0.7625685	0.465
Residual	47	1.3271997	0.924963	NA	NA
Total	52	1.4348679	1.000000	NA	NA

pairwise<-pairwise.adonis(beta_div_func_post5$S,post5_arrange$type_time, perm=999)
pairwise%>%
  tt()

pairs	Df	SumsOfSqs	F.Model	R2	p.value	p.adjusted	sig
Control.5_Post-FMT1 vs Treatment.5_Post-FMT1	1	0.1195408549	4.84764704	0.2442429086	0.067	1
Control.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0525878365	1.77308932	0.0997625840	0.264	1
Control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0265995825	1.17541806	0.0684360667	0.312	1
Control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0145818992	0.69975992	0.0419023938	0.388	1
Control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0080695208	-0.21617323	-0.0136958691	0.923	1
Treatment.5_Post-FMT1 vs Hot_control.5_Post-FMT1	1	0.0129803540	0.44307662	0.0286909552	0.456	1
Treatment.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	0.0267162134	1.22560581	0.0755352882	0.337	1
Treatment.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0384388433	1.93281582	0.1141461550	0.218	1
Treatment.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	0.0553988290	1.47819391	0.0897060633	0.250	1
Hot_control.5_Post-FMT1 vs Treatment.6_Post-FMT2	1	-0.0040061386	-0.14850469	-0.0093684974	0.729	1
Hot_control.5_Post-FMT1 vs Control.6_Post-FMT2	1	0.0024023972	0.09538980	0.0059265296	0.612	1
Hot_control.5_Post-FMT1 vs Hot_control.6_Post-FMT2	1	-0.0004960759	-0.01190328	-0.0007445087	0.842	1
Treatment.6_Post-FMT2 vs Control.6_Post-FMT2	1	-0.0080428882	-0.44298625	-0.0284750185	0.838	1
Treatment.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	-0.0011796256	-0.03404738	-0.0021324990	0.912	1
Control.6_Post-FMT2 vs Hot_control.6_Post-FMT2	1	0.0036300838	0.11048757	0.0068581148	0.681	1

beta_richness_nmds_post5 <- beta_div_richness_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_neutral_nmds_post5 <- beta_div_neutral_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_phylogenetic_nmds_post5 <- beta_div_phylo_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

beta_functional_nmds_post5 <- beta_div_func_post5$S %>%
                metaMDS(.,trymax = 500, k=2, verbosity=FALSE) %>%
                scores() %>%
                as_tibble(., rownames = "sample") %>%
                left_join(post5_nmds, by = join_by(sample == Tube_code))

p0<-beta_richness_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Richness beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p1<-beta_neutral_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
        scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y = "NMDS2", x="NMDS1 \n Neutral beta diversity") +
                theme_classic() +
                theme(legend.position="none")
  
p2<-beta_phylogenetic_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Phylogenetic beta diversity") +
                theme_classic() +
                theme(legend.position="none")

p3<-beta_functional_nmds_post5 %>%
            group_by(type, time_point) %>%
            mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
            mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
            ungroup() %>%
            ggplot(., aes(x=NMDS1,y=NMDS2, color=type, shape=time_point)) +
                scale_color_manual(name="Type",
                       breaks=c("Control", "Hot_control", "Treatment"),
                       labels=c("Cold-Cold", "Hot-Hot", "Cold-Hot"),
                       values=c("#4477AA","#d57d2c","#76b183")) +
                geom_point(size=2) +
                geom_segment(aes(x=x_cen, y=y_cen, xend=NMDS1, yend=NMDS2), alpha=0.2) +
        labs(y= element_blank (), x="NMDS1 \n Functional beta diversity") +
                theme_classic()+
                theme(legend.position="none")

ggarrange(p0, p1, p2, p3, ncol=2, nrow=2, common.legend = TRUE, legend="right")